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Study of the Hubbard model on the triangular lattice using dynamical cluster 
approximation and dual fermion methods 
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We investigate the Hubbard model on the triangular lattice at half-filling using the dynamical 
cluster approximation (DCA) and dual fermion (DF) methods in combination with continuous-time 
quantum Monte carlo (CT QMC) and semiclassical approximation (SCA) methods. We study the 
one-particle properties and nearest-neighbor spin correlations using the DCA method. We calculate 
the spectral functions using the CT QMC and SCA methods. The spectral function in the SCA and 
obtained by analytic continuation of the Fade approximation in CT QMC are in good agreement. 
We determine the metal-insulator transition (MIT) and the hysteresis associated with a first-order 
transition in the double occupancy and nearest-neighbor spin correlation functions as a function 
of temperature. As a further check, we employ the DF method and discuss the advantages and 
limitation of the dynamical mean field theory (DMFT), DCA and recently developed DF methods 
by comparing Green's functions. We find an enhancement of antiferromagnetic (AF) correlations 
and provide evidence for magnetically ordered phases by calculating the spin susceptibility. 

PACS numbers: 71.10.Fd 
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I. INTRODUCTION 

The physics of systems which exhibit strong electronic 
correlations and geometric frustration at the same time 
is still unclear and therefore interesting. Recent exper- 
iments, such as discovery of the pyrochlore compound 
LiV204^ which show heavy fermion behavior and organic 
materials k-(BEDT-TTF)2X- which exhibit various in- 
teresting phases, motivated us to study the frustrated 
systems in more detail. Theoretically they are described 
by a two-dimensional one-band Hubbard or t-J models 
on the triangular lattice. It is well known that on the 
square lattice at half-filling the ground state is a Mott 
insulator with AF order but on the triangular lattice the 
frustration suppresses AF order and we expect to find a 
Mott transition. 

In this paper, we study the model which was presented 
in a recent paper of Imai and Kawakami^. They used 
the DCA method^^'^® in combination with noncrossing 
approximation (NCA) and fluctuation exchange (FLEX) 
methods at high temperature regions in metallic states to 
demonstrate how geometrical frustration suppresses AF 
correlations by tuning aniotropic hopping t' in Fig. 1(a). 
However, the methods used in the paper are limited to 
high temperatures. For this reason we investigate the 
low temperature MIT with a first-order transition and 
the evidence of magnetically ordered phases. In addition, 
we test the newly developed DF metho d^'^'^^i^^ beyond 
the single-site DMFT method by comparing the Green's 
function of single-site DMFTii^, DF and DCA methods 
with Nc=A and Nc=l&. 

The paper is organized as follows: In Sec. [TTl we in- 
troduce the model and discuss the advantages and lim- 
itations of the computational methods briefly. In Sec. 
IIII[ we present the results. In the first part, we com- 
pare the spectral functions which exhibit the quasipar- 
ticle peak and gap structure obtained by the SC A^'^i^-'^ 



and CT QMG^^ methods with Fade approximation. In 
the second part, we show that the MIT is a first-order 
transition by measuring the total density of state (DOS), 
double occupancy and nearest-neighbor spin correlations. 
In the third part, we calculate the single-particle Green's 
function using the DMFT, DF and DCA methods with 
Nc—'^, and Nc—1&- Specifically, we show that the DF 
method which is based on the single-site DMFT method 
can describe non-local correlation effects very well. In 
the last part, we explore the spin susceptibility using the 
DF method. Sec llVi we give a summary work. 



II. MODEL AND METHODS 

A. Model 

We consider the two-dimensional Hubbard model on 
the triangular lattice. 



^ = -iE^U- + f^E 
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where Cio-(cJ^) is the annihilation (creation) operator of 
an electron with spin a at the i-th site, t is the hopping 
matrix element and U represents the Coulomb repulsion. 
In this paper we only consider the isotropic hopping of 
t'=t in Fig. l(a)-(b). 

Due to the geometrical frustration, this model has bro- 
ken particle-hole symmetry even at half-filling unlike the 
case of the square lattice and the original Brillouin zone 
(BZ) has an hexagonal structure shown by the doted line 
in Fig. 1(c). For this model there are a lot of stud- 
ies using a variety of methods such as the path inte- 
gral renormalization group method^, the quantum Monte 
carlo method^, the DMFT^i^ and the cluster extension 
method of DMFT^iLlLi^J^. Especially, DFiiii^i^ and 
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FIG. 1: (a) Schematic representation of triangular lattice with 
electron hoppings. (b) Equivalent representation of (a) for a 
square structure, (c) Example of the coarse-graining cells in 
the BZ for the triangular lattice (a), where the cluster size is 

Nc = 4. 



DCA methodsiiii^ are noteworthy because both meth- 
ods can capture non-local correlation effect which are 
lost in the single-site DMFT and they are computational 
cheaper and have less of a sign problem than the lat- 
tice QMC calculation. In short, DCA method can treat 
correlations up to a cluster size Nc accurately and long 
range correlations are considered on the mean- field level. 
On the other hand, the DF method considers long range 
as well as short range correlations within pertubative di- 
agram expansion, which is done by introducing an aux- 
iliary field. Because each method has its limitations, it 
is useful to compare results of both. We use the CT 
QMC methodic, which can access the low temperature 
region easily without the Trotter error, as well as the 
SCA method^Si^ as impurity solvers. 



B. DCA method 

The DCA metho d^^'^^ assumes that the self-energy in 
the first BZ is constant and the coarse-grained Green's 
function (DCA equation) is given by Eq. (2). 



G.(K,z) 



-E- 



K 



^K+K 



S.(K,z)' 



(2) 



where N is the number of lattice sites in each first BZ 
and the summation over K is calculated in each of them. 
For delimitation we consider an example of Nc—4: in or- 
der to explain the DCA method. The first BZ (dashed 
line in Fig. 1(c)) is created by partitioning the original 
BZ. Like the standard procedure of DMFT, the coarse- 
grained Green's function is determined self-consistently 
after several iterations. The main advantages of the DCA 
method are that it considers short range correlations in 
the cluster size exactly and has smaller computational 



load and fermionic sign problem compared to the lat- 
tice calculation by QMC method. However, it is still 
expensive in terms of computational time and long range 
correlations are just treated on the mean-field level. 



C. DF method 

The DF methodi^ii^ii^ is a relatively new method 
which can describe non-local correlations based on the 
single-site DMFT method. The basic idea of the DF 
method is to convert the hopping of different fermions 
into an effective coupling to an auxiliary field. Each lat- 
tice site can be viewed as an impurity which is easily 
described by the DMFT method. While these impurities 
are not totally isolated, they are perturbatively coupled 
by the auxiliary field. The starting point is the action of 
DMFT which is represented in the form 



^[c+,c]=^SLp-E(^''- 



(3) 



where Aj/ is the hybridization function describing the in- 
teraction of an effective impurity with a bath and u is the 
fermionic Matsubara frequency. Here we use the disper- 
sion relation Ck = —2t[cos{kx) + cos{ky) + cos{kx + ky)] 
based on the correspondence of a triangular lattice to a 
square lattice with diagonal hopping (b) in Fig. 1. This 
"square lattice" has a simple BZ which makes the mo- 
mentum summation to be easily performed by using the 
fast Fourier transformation (FFT). By the dual transfo- 
mation, the lattice problem is changed to an impurity 
problem which is coupled by the auxiliary field /(/^). 
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The lattice Green's function is derived from the exact 
relation between Eq. (3) and Eq. (4). 
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where Gf, j. is the dual Green's function and g^ is the 
local Green's function calculated by single-site DMFT. 
The main point for this method is that in Eq. (4) the 
integration over c^ and c can be performed separately for 
each site which yields an effective action of the auxiliary 
field / and p . The Taylor expansion in powers of /^ 
and / will introduce the two, three, . . . ,-particle ver- 
tex functions. Using the skeleton-diagram expansion we 
calculate the dual self-energy and dual Green's function 
by the Dyson equation. We obtain the lattice Green's 
function via Eq. (5). Even though the DF method is 
an approximate method, it considers not only the short 
range but also the long range correlations. Moreover, the 



calculation of the two-particle properties does not intro- 
duce serious computational burden and fermionic sign 
problem. 



D. SCA method 

At high temperature the Monte carlo integration over 
the auxiliary classical field (J>{t) can be approximated by 
assuming 4>{t) « const. This approximation is useful be- 
cause it allows to check the QMC results at temperature 
quickly. In this part we introduce the SCA metho d^"'^^ 
as impruity solver for DCA method. We consider a four- 
site cluster (7Vc=4) for triangular lattice like the structure 
of Fig. 1 (a) . In this case the partition function is defined 
as a functional integral over 2 x 4-component spin and 
site-dependent spinor fields c^ and c as 



where M is defined as 
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Here a.cr{T,T') is the Weiss field which is determined 
self-consistently by Eq. (2) and /3=1/T is the inverse 
temperature. In this model a.^ is given as 



aaiT,T') = 
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We can decouple the interaction term as 

C/n,t(r)n,i(T) = ^ [iVf (r) - Mfir)] , 



(8) 



with n^ni = \ ((n-r + n^f - (n^ - n^f) = \{N''~M'^). 
We employ the continuous Hubbard-Stratonovich(IIS) 
transformation in order to decouple M terms related to 
auxiliary field 4>j(t\ Here we assume that (?!>j(t) is r in- 
dependent and N term is neglected because charge fiuc- 
tuations are small at half-filling. By a Grassmann in- 
tegration we can rewrite the partition function which is 
represented as a four-dimensional integration in terms of 
(\)j and the fermionic Matsubara frequency. 
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where the effective action 5*6// = PV is defined by 
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Here j=l,2,3,4 and cr^ is the z-component Pauli matrix. 
The impurity Green's function is calculated by 



G,= 
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In real frequency space the spectral function is calculated 
by replacing a(iw) to a(w). The SCA method is not only 
cheap in computational time but also gives good results 
in the strong coupling regime. On the other hand, it 
underestimates the spectral function around lj — Q and 
gives qualitatively wrong results at low temperature. 



E. CT QMC method 

Here we describe the CT QMC method^^. The starting 
point is action can be split into an unperturbed action 
S{) and an interaction part W. By Taylor expansion of 
partition function in powers of the interaction U, we can 
reexpress the partition function 
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where Zq ~ Tr{Te^^°) is the partition function for un- 
perturbed system, integration over dr implies the inte- 
gral over T and sum over all lattices states and T is the 
time-ordering operator. By Wick's theorem the weight 

is determined by 



function Drl'.'.'.rk 
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,) is the bare Green's function. The 



Green's function is defined by 
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and by the fast-update formul a^^i^^ and the Fourier 
transformation we can rewrite the Green's function in 
the Matsubara frequencies space. 

G{u:)=go{u^)~go{uj)[-Y,M.,,e^-^-^-^^^]g^{u:) (17) 

where M — D~^ and go{uj) is the bare Green's function. 
A two-particle Green's function related to a vertex func- 
tion for DF method can be calculated by Wick's theorem. 
With this method it is possible to perform calculations in 
low temperature regions which cannot be accessed easily 



with determinant QMC method without Trotter error. 
For example, the matrix size of the CT QMC method is 
scaled by k ~ 0.5NJJf3 which is comparable to determi- 
nant QMG22. method scaled by k ~ bNJjp. Moreover, 
even if the recently developed strong-coupling CT QMC 
method )^^'^^ which is based on a diagrammatic expan- 
sion in the impurity-bath hybridization, has nice advan- 
tages such as removing the fermionic sign problem and 
calculating lower temperature regions, its computational 
effort in large cluster is increased exponentially by the 
number of sites Nf,- However, our CT QMC method can 
overcome the problem because the computational burden 
only increases linealy with the number of sites. 



III. RESULTS 

A. Comparison of the spectral functions for the 
SCA and C-T QMC methods 

First we compare the one-particle spectral functions 
obtained from the SCA and CT QMC methods with Fade 
approximation for analytical continuation. Since the pro- 
cess, in which G{iuj) calculated by the CT QMC method 
changes into G{uj) with Fade approximation, introduces 
large error, it is useful check to compare QMC results to 
SCA results which are calculated in real frequency space. 
Moreover, because the systems with geometrical frustra- 
tion have large Uc, the SCA method is suitable for this 
model. The spectral function is given by 



A,{K,u;) = --lmG,iK,cj) 



(18) 



We compare the Green's function on the Fermi surface 
K=(7r, 7r/-\/3). The results are shown in Fig. 2(a)-(b). 
At U=6 in Fig. 2(a) the difference of both results is that 
the peak of quasiparticle obtained from QMC method 
lies around the Fermi level(Lj=0) due to the geometrical 
frustration. On the other hand, the peak obtained from 
SCA method deviates from the Fermi level because the 
SCA underestimates the a;=0 peak. At U=9 in Fig. 2(b) 
the agreement of both results is more reasonable and a 
(pseudo)gap structure is represented. 
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FIG. 2: One-particle spectral function A(K, to) corresponding 
to the K=(7r,7r/y3) for /3 = 1.6667, (a)U=6 and (b)U=9 by 
means of the SCA and CT QMC with Pade approximation. 



is a small gap at A'c=16 which completely disappears 
at A^c=64. In this system we did not find a band insu- 
lator transition. However, on the triangular lattice the 
frustration is enough to destroy the AF correlation. In 
Fig. 3(a)-(b), we can see the MIT by comparing the total 
DOS for U=6, U=IO and /3=4 using the DCA method 
with iVc=4 and Nc=l&- Unlike the results for unfrus- 
trated square lattice, the quasiparticle peak around the 
Fermi level is clearly seen with increasing Nc at U=6 in 
Fig. 3(a). At U=10 in Fig. 3(b) we can see the Mott 
insulator in both Nc—^ and Nc—IQ- This is evidence of 
a MIT on the triangular lattice. In the low temperature 
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B. The metal- insulator transition with a first-order 
transition 



FIG. 3: Total DOS with N^^A and Nc^l& for /?=4, (a)U=6 
and (b)U=10 via CT QMC with Pade approximation. 



Here we present our results on the MIT due to geo- 
metrical frustration effect obtained with the CT QMC 
method. In previous study of unfrustrated square lattice 
using DCA method, it was shown that short-range AF 
correlations destroy the Fermi liquid quasiparticle peak 
at finite temperature^^. According to Ref. 25, the au- 
thors increased the system size gradually in the weak- 
coupling regime on unfrustrated square lattice and mea- 
sured the total DOS. Eventually, even if the quasiparticle 
peak is clearly visible at Nc=l in DMFT method, there 



regime we are also interested in finding whether there is a 
first-order transition or a continuous transition and how 
the geometrical frustration effects the system. We expect 
our system to have a first-order transition because of the 
recent two cellular DMFT (CDMFT) resultsi^iSe which 
show a first-order transition on the square lattice with 
Nc—'^ and on the triangular lattice with anisotropic hop- 
ing at low temperatures. In order to find evidence of a 
first-order transition we measure the double occupancy at 
several temperatures. Our result is shown in Fig 4. The 



system displays a crossover between metal and insulator 
at T=0.2. At T=0.1 we can see hysteresis associated with 
a first-order transition and at lower temperature hystere- 
sis is more clear. In order to understand the system more 
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FIG. 4: Double occupancy as a function of U/t at several tem- 
perature. )7c=7.2 for T=0.2, [/c=6.9 for T=0.1 and Uc=6.7 
for T=0.05. 

clearly we calculate the nearest-neighbor spin correlation 
function {S^Sfj^^ which is shown in Fig. 5. At Uc jumps 
of the spin correlation function indicate the MIT aris- 
ing from competetion between the quasiparticle forma- 
tion and the frustrated spin correlation. Specifically, the 
spin correlation is enhanced weakly at Uc— 7.2 for T=0.2 
while it is increased rapidly at Uc in T < 0.1. Here is 
Uc^6.9 for T==0.1 and C/c=6.7 for T=0.05. This means 
that the entropy at T=0.2 and in T < 0.1 is released by 
geometrical frustration and spin correlation, respectively 
as temperature decreases and the entropy at insulator 
state in T < 0.1 has small value which is triggered a 
first-order transition because of a formation of AF state. 
Moreover, we find that the anomalous character in the 



state the spin correlation is weak. This is the reason that 
geometrical frustration is more dominant than AF spin 
correlation at lower temperature in the metallic state be- 
cause the frustration on the triangular lattice is stronger 
than that on the Kagome lattice. However, in the in- 
sulating state AF spin correlation is enhanced stronger 
than the frustration effect at lower temperature so the 
spin correlation is strong with decreasing temperature. 



C. Comparison of Green's functions among the 
DCA, DF and DMFT methods 

In this part wc using the DMFT, DF and DCA meth- 
ods with Nc—4: and Nc=16 to study the non-local correla- 
tion effects and compare the on-site and nearest-neighbor 
Green's functions in the Matsubara space. 




FIG. 6: The imaginary part of on-site Green's function for 
/3=4, (a) U=6 and (b) U=10. The real part of nearest- 
neighbor Green's function for /3=4, (c) U=6 and (d) U=10. 
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FIG. 5: The nearest neighbor spin correlation function as a 
function of U/t at several temperature. Uc=7.2 for T=0.2, 
Uc^6.9 for T=0.1 and [7^=6.7 for T=0.05 

metallic state is unlike the results of the nearest-neighbor 
spin correlation on the Kagome lattice^^. In the metallic 



In Fig. 6(a)-(d), we present the Green's functions ob- 
tained from DMFT, DF and DCA method with Nc^i 
and Nc=16 for /3=4, U=6 and U=10. The on-site Green's 
function of DMFT method in Fig. 6(a) is similiar to the 
results of DCA and DF method and all of these indi- 
cate the metallic states. A remarkable point is that in 
Fig. 6(a) and (b), both the on-site and nearest-neighbor 
Green's functions obtained from DF method are closer 
to those of the DCA method with 7Vc=16 than Nc='i. In 
Fig. 6(c) at U=10 the on-site Green's function calculated 
by the DMFT still shows the metallic state which over- 
estimates the value of Uc because of a lack of non-local 
correlation. However, the DCA and DF methods can 
capture the insulating state and the agreement of the 
on-site Green's function calculated by the DF and DCA 
methods with 7Vc=16 is quite reasonable. In Fig. 6(d), 
the nearest-neighbor Green's functions obtained from DF 
method are still closer to those of Nc=16 than Nc=4:. 
This suggests that despite the fact that the DF method 
is a perturbative method, it would describe physics quite 
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FIG. 7; (a) The spin susceptibility xil) in the insulating state 
for U/t=10.0 and f3t = 2.5. (b) The spin susceptibility as a 
function of temperature at q=(0,0) and q = {2n/3, 2tt/3). 



well than the DCA method with small cluster size. We 
expect that considering high-order diagrams will improve 
the results of the DF method. 



D. The spin susceptibility using the DF method 

In order to explore a magnetic instability we measure 
the spin susceptibility using the DF method. The reason 
why we employ DF method for the spin susceptibility is 
that the cluster-extension method of the DMFT takes a 
large amount of time in order to obtain the two-particle 
properties. On the other hand, because the DF method 
includes the vertex renormalization through the Bethe- 
Salpeter equation, the computational burden is not seri- 
ous and the results are relatively good compared to those 



of QMC methodic. Fig. 7(a) shows x(g) for U/t=10.0 
and (3t — 2.5 where the system is in the insulating state. 
The x{q) has a maximum peak at q = (27r/3, 2tt/3). The 
spin susceptibility xil) a.t q — (27r/3, 27r/3) and q=(0,0) 
as a function of temperature is exhibited in Fig. 7(b). As 
temperature decreases, xil) 3± q = (27r/3, 27r/3) shows 
strong enhancement of the AF correlations. 



IV. CONCLUSIONS 

In summary, we have investigated the Hubbard model 
on the triangular lattice using the DCA and DF method. 
Using the DCA method we compared the spectral func- 
tions obtained from SCA and CT QMC methods. We 
found a good agreement of both methods and the quasi- 
particle peak and gap structure are presented in the weak 
and the strong coupling regions, respectively. We found 
a MIT with a first-order transition at low temperatures 
because of the effect of geometrical frustration. More- 
over, we employed the DF method which considers the 
long range as well as short range correlations and com- 
pared the Green's functions of the DF method to those 
of the DMFT and DCA method with Nc^i and A^c=16. 
We found that the DF method does not only overcome 
the overestimation of Uc in DMFT method but also that 
its results are closer the case of Nc=16 than A^c=4. Fi- 
nally, we calculated the spin susceptibility xil) via DF 
method. We found that the xil) at g = (27r/3, 27r/3) 
grows rapidly as temperature decreases. 

We would like to thank A. Lichtenstein and E. Gorelov 
for their assistance in implementing the CT QMC code. 
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